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A  HYBRID  METHOD  FOR  PARAXIAL  BEAM 


PROPAGATION  IN  MULTIMODE  OPTICAL  WAVEGUIDES 

by 

David  C.  Chang  and  Edward  F.  Kuester 

Electromagnetics  Laboratory 
Department  of  Electrical  Engineering 
University  of  Colorado 
Boulder,  Colorado  80309 

Abstract 

A  hybnid  [non-nay,  non-modat )  method,  ion  Computing  the  iieldt 
oi  a  panaxial  beam  pnopagating  in  a  multimode  waveguide  (panallei- 
plate  on  dielectnic  tlab]  at  Zange  axial  distances  it  pnetented. 
The  method  ii>  bated  on  the  Founien  and  F netnei  teli~  imaging 
pnopentiet  o&  thete  waveguide* ,  and  it  capable  oi  high  accunacy. 
The  method  it  much  mone  evident  than  nay  on  mode  appnoachet, 
while  giving  complete  iield  inionmation  which  coupled-powen 
equation t  do  not  pnovide. 


I.  INTRODUCTION 


Multimode  optical  fibers  appear  at  present  to  be  the  most  common 
optical  waveguiding  medium  for  applications  in  the  immediate  future. 
Incoherent  sources  and  relatively  simple  detectors  can  be  used,  and 
the  tolerance  problems  encountered  with  single-mode  fibers  are  far  less 
severe  with  such  waveguides. 

At  present,  there  are  essentially  three  methods  available  for  field 
computation  in  multimode  waveguides.  First,  one  can  take  a  pure  modal 
approach—the  excitation  amplitude  of  each  mode  is  computed,  anu  all 
modes  are  summed  together.  Although  in  principle  exact,  this  approach 
suffers  not  only  from  the  large  number  of  modes  which  must  be  kept  track 
of  ( 1 00 ^ 1 000  for  a  typical  fiber;  30 ^ 1 00  for  a  slab  geometry)  but  also 
from  a  large  degree  of  cancellation  of  terms  in  the  mode  sum  when  the 
field  does  not  match  that  of  an  individual  mode.  Examples  of  the  appli¬ 
cation  of  this  method  may  be  found  in  [1].  Although  in  some  special 
cases  approximate  closed-form  results  are  available,  a  computer  analysis 
is  generally  required,  and  roundoff  errors  can  be  expected  to  accumulate, 
especially  for  large  propagation  distances. 

A  second  approach  is  that  of  geometrical  optics  (sometimes  encountered 
as  the  WKB  method).  An  excellent  discussion  of  this  approach  has  been 
given  by  Gloge  and  Marcatili  [2]  (see  also  [3]).  Here  one  approximates 
the  effect  of  a  large  number  of  discrete  propagating  modes  by  a  continu¬ 
ously  distributed  propagation  constant  belonging  to  a  "continuous  spectrum" 
of  modes.  These,  when  computed  under  the  WKB  approximation,  can  be 
interpreted  as  a  cone  of  rays  lying  within  some  characteristic  acceptance 
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angle  of  the  fiber.  The  propagation  problem  then  reduces  to  that  of 
determining  the  amplitude  with  which  each  >'ay  is  excited,  and  tracing  it 
down  the  length  of  the  guide.  Intuitively  more  suitable  for  multimode 
guides  because  of  the  "high-frequency"  nature  of  the  problem,  this  approach 
is  nonetheless  approximate  by  virtue  of  the  geometrical  optics  technique. 
Moreover,  in  a  situation  where  paraxial  propagation  conditions  exist  (see 
below),  a  large  number  of  rays  can  be  expected  to  contribute  at  large 
propagation  distances  (hundreds  of  meters  or  several  kilometers  may  not 
be  uncommon).  In  this  region,  the  geometrical  optics  approach  can  be 
seen  to  suffer  from  similar  disadvantages  as  does  the  first. 

A  third  approach  (see,  e,g.,  [4], [5])  is  a  purely  numerical  one, 
wherein  the  partial  differential  equation--Helmholtz  or  its  parabolic 
approximant— is  tackled  directly,  without  the  use  of  either  mode  or  ray 
concepts.  In  [4]  and  [5],  the  equation  is  discretized  and  solved  with 
the  aid  of  fast  Fourier  transform  techniques.  This  method,  like  the  first, 
is  also  capable  of  arbitrary  accuracy  in  principle,  and  requires  neither 
a  detailed  knowledge  of  a  large  number  of  modes,  nor  the  tracing  of  a 
large  number  of  ray  paths.  Again,  however,  when  very  long  propagation 
distances  are  being  studied,  the  discretization  of  the  wave  equation  in 
the  longitudinal  direction  can  lead  to  large  error  accumulations  which  do 
not  seem  easily  avoidable  by  this  technique. 

Finally,  we  might  also  mention  here  the  coupled-power  equations 
approach  [1].  This  method  seeks  only  to  find  the  total  power  carried  by 
each  mode,  since  for  many  applications  the  details  of  the  field  distribu¬ 
tion  from  each  mode  are  not  of  interest.  One  then  takes  a  statistical 
approach  to  these  equations,  and  obtains  useful  results  for  pulse  dispersion 


when  each  mode  of  the  guide  is  detectable  only  through  its  total  power. 

There  are  many  other  applications,  however,  when  the  fields  themselves  are 
important,  such  as  in  the  design  of  couplers,  splitters,  switches,  splicers, 
etc.,  and  it  is  this  problem  in  which  we  are  interested  here. 

The  method  we  propose  is  based  on  the  imaging  properties  of  multi- 
mode  waveguides.  In  the  paraxial  approximation,  a  parallel-plate  or 
dielectric  slab  waveguide  will  periodically  reconstruct  the  field  pattern 
at  the  input  plane  (and,  at  more  frequent  intervals,  a  string  of  such 
replicas).  Because  of  this,  we  need  only  perform  our  field  computations 
within  the  space  of  one  of  these  periods,  and  will  not  suffer  the  loss  of 
accuracy  at  large  distances  associated  with  the  methods  described  above. 

Our  computations  will  be  performed  for  a  parallel -pi ate  waveguide  with 
perfectly  conducting  walls,  but  the  results  are  immediately  applicable  to 
the  dielectric  slab  waveguide  (Appendix  A).  The  method  will  allow  a 
simple  formula  to  be  obtained  for  the  propagation  of  a  Gaussian  beam  of 
substantially  narrower  width  than  that  of  the  guide.  (Note  that  the  study 
by  Felsen  and  Shin  [6]  of  beam  propagation  in  waveguides  is  in  practice 
restricted  to  beams  which  propagate  obliquely  to  the  axis  of  the  guide, 
and  actually  will  suffer  from  the  same  drawbacks  at  large  axial  distances 
as  does  the  ray  method,  although  it  is  probably  superior  to  the  ray  method 
at  shorter  distances).  Numerical  comparisons  with  exact  (modal)  calcula¬ 
tions  are  quite  favorable,  even  for  rather  large  propagation  distances. 
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II.  THE  PARAXIAL  APPROXIMATION 

To  fix  ideas,  let  us  consider  the  parallel-plate  waveguide  illustrated 

in  Fig.  1  (the  discussion  of  this  section,  however,  is  quite  general  and 

need  not  be  restricted  to  this  specific  waveguide).  The  walls  at  x  =  0 

and  x  =  a  are  perfectly  conducting,  and  some  known  source  produces  a  given 

excitation  or  input  field  at  the  plane  z  = 0.  For  simplicity,  we  restrict 

ourselves  to  two  dimensional,  TE  fields,  so  that  the  entire  field 

H  =  a  H  +  a  H  ,  E  =  a  E  ,  where  a  ,a  ,5  are  Cartesian  unit  vectors,  can 
xx  z  z  y  y  xyz 

be  derived  from  the  scalar  function  E^  which  satisfies 

{-4*  "4  +  k2)Ey  *  0  (,) 

'  8x2  dz  y 

for  z>0.  Here  k  =  w/fte,  where  a  time  dependence  exp(iwt)  has  been 
assumed,  and  y,  e  are  the  electrical  parameters  of  the  medium  filling 
the  waveguide. 

In  the  paraxial  approximation,  we  write 

Ey( x, z)  =  e"lkzA(x,z)  (2) 

and  assume  that  most  propagation  takes  place  nearly  in  the  z-direction; 
that  is,  A(x,z)  as  a  function  of  z  varies  slowly  compared  to  exp(-ikz). 
Inserting  (2)  into  (1)  we  obtain 

( “  r  -  2ik  )  A(x,z)  =  0  (3) 

3XZ  92  9z 

2  2 

and  in  the  paraxial  approximation,  we  neglect  the  9  /9z  term  compared 
to  the  first  derivative  term  because  of  the  slow  variation  of  A(x,z)  in 
z  assumed  above.  We  thus  obtain  the  following  parabolic  equation  for 
A(x,z)[7]: 
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(l7~ 21k  3z 


)  A(x,z)  = 


0 


(4) 


To  put  this  approximation  on  a  more  quantitative  footing,  we  can  apply 
some  ideas  from  the  boundary-layer  technique  [8].  By  "stretching"  the 
variable  x  into  a  new  variable  v  =  x/a,  we  can  deal  with  a  transverse 
variable  v  which  is  0(1)  over  the  entire  cross-section  of  the  guide.  When 
rewritten  in  terms  of  v,  equation  (4)  becomes 

ii 


|  -  2ika2  ~  j  A(v,z)  =  0 


(5) 


which  suggests  that  the  scaling  c  =  const  z/ka  might  be  convenient.  For 

7 

reasons  which  will  become  clear  in  the  next  section,  we  choose  ?  *  nz/4ka . 
Making  both  changes  of  variable  in  (3),  we  have 

flL.  4*  JL.Wd  -  -  — \ 

W  T*>  ^ 

It  is  natural  now  to  assume  a  solution  to  (6)  of  the  form 


(6) 


A(v,s)  -  A 0(v,s)  +  A j(v,c)  +  -5-3-  A2(v,c)  + 
k  a  k  a 

_2 

so  that,  by  matching  powers  of  (ka)  ,  we  find  a  recurrent  set  of 
equations  for  AQ,  A^ ,  etc.: 


(7) 


(^'T  w)Ao('',s)  "  0 

2  \  ^ 

[  -3*-  _  TIL  _S.  )  A,(v»c)  '  “ 

l  J  2  '  16 


2  9  Ao(v’  d 


(8) 

(9) 


and  so  on.  As  initial  conditions  at  t,  =  0,  we  require  AQ(x/a,0)  =  E^(x,0)- 


a  given  function— while  A^U.O)  =  A^W.O) 


...  *  0. 
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We  are  now  in  a  position  to  estimate  the  magnitude  of  the  correction 
term  A^,  and  therefore  the  error  involved  in  the  paraxial  approximation, 
tt  is  easily  verified  that,  for  example, 


A ,,  ;1  ■  is  3V-S)  -  i.c 

VV,U  2tt  ,..4  8  7JT 


(10) 


3v  w  dC 

with  A^,  Aj,  etc.,  given  by  similar  expressions.  Because  of  the  scaling 

of  the  variables  v  and  c,  the  differentiations  in  (10)  do  not  increase 


thi>  order  of  magnitude  of  the  function  Aq.  Thus  if  Aq  is  0(1),  then 


I A-]  (v,s)  |  <  0(5) 


(ID 


so  a  :riterion  for  the  accuracy  of  the  paraxial  approximation  (A  -  A  ) 

2  2 

is  that  c/k  a  «  1,  or  in  other  words, 

kz  «  (ka)4  (12) 

2  2 

In  addition,  of  course,  we  also  have  the  condition  k  a  »  1,  which  is 
implicit  in  the  expansion  (7)  and  the  fact  that  the  guide  is  highly  multi- 
mode. 

These  arguments  are  not  restricted  to  the  case  of  a  parallel -pi ate 
waveguide.  However,  as  with  any  "order-of-magnitude"  arguments,  they 
say  nothing  about  the  proportionality  constant  implicit  in  (11).  In  fact, 
this  constant  will  depend  sensitively  on  the  function  Aq(v,c),  mostly 
through  its  initial  value  Aq(v,0).  Thus,  a  detailed  study  of  the  esti¬ 
mate  (11)  should  be  made  when  using  the  paraxial  approximation  in  any 
specific  situation.  This  procedure  is  discussed  in  some  detail  by 
Tappert  [9],  who  also  gives  a  large  number  of  references  to  its  use  in 
acoustics.  In  Appendix  A,  we  consider  the  paraxial  approximation  for 


a  dielectric  slab  waveguide,  and  demonstrate  an  approximate  equivalence 
wtth  a  parallel -pi ate  waveguide. 

In  closing,  we  might  also  note  that  Polyanskii  [11]  has  obtained  a 
formula  relating  the  solution  of  the  parabolic  equation  to  that  of  the 
Helmholtz  equation  which  is  an  alternative  to  the  perturbation  series 
(7).  Further,  another  variant  of  the  parabolic  equation  more  suitable 
for  off-axis  propagation  has  been  proposed  in  [12]. 
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III.  GREENES  FUNCTION  AND  IMAGING 


By  well-known  techniques,  the  field  E^(x,z)  for  z  >0  in  the  waveguide 
can  be  expressed  in  terms  of  the  field  Ey(x,0)  at  an  input  plane  (z  =  0) 


by  means  of  a  Green's  function  G(x,x‘;z): 


a 

E  (x,z)  =  E  (x* ,0)G(x,x‘ ;z)dx' 

y  Jn  y 


where  G  can  be  expressed  as  a  modal  expansion: 

G  (x,x';2)-i  J  sin  »  sin  ’sf- z  >  0  (14) 

3  q  a  a  — 

m- 1 

where  =  (k  *mS  /aj  .  On  the  other  hand,  the  Green's  function  Gq 
for  the  paraxial  approximation  (2),  (4)  to  Ey  is  given  by: 

G  (x.x'u)  =  ^e-,k2  l  sjnap.s1nSIiIeiz"2,'2/2ka2  (15) 

o  a  a  a 

m- 1 

Evidently  Gq  could  have  been  obtained  from  G  by  replacing  by  the 

2  2  2 

first  two  terms  of  its  binomial  expansion  Gm  =  k-m  tt  /2ka  .  We  see  that 
whereas  G  has  a  large  but  finite  number  of  modes  which  propagate,  the 
paraxial  approximation  Gq  has  infinitely  many  such  modes.  It  is  con¬ 
venient  to  rewrite  Gq  by  expressing  the  sine  functions  as  exponentials; 
the  result  is: 

1 

G  (x.x1  ;z)  =— e"ikz  £  eizmir/2ka ,  e-iimr(x-x'  )/a  _  e-iimr(x+x'  )/a>  (16) 

^  ^  m=  -oo  J 


If  we  define 


z^  =  4ka  /it 


(so  that  the  stretched  variable  of  the  previous  section  is  5  =  z/z^),  then 
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.  i  ®  2wt«fz/z11 

Go(x*x  ;z)  =  2?  e  ik  i  e  11 < 


e-iiwr(x-x‘)/a  _e-inm(x+x')/al  ^18j 


Equation  08}  is  the  basis  for  the  so-called  Fourier-and  Fresnel- 
imaging  properties  of  this  waveguide  [13],  It  is  easily  seen  from  (18) 
that  exp(ikz)G0  is  a  periodic  function  of  z: 

SgCx.x' ;z  +  )  *  G0(x,x';z)e  11  (19) 


In  particular,  since  (13)  implies  that  G(x,x*;0)  is  equal  to  fi(x-x') 
for  0<(x,x')<a,  we  have 


Go(x,x';  nz^) 


-iknz 

S(x-x*)e 


11 


(20) 


for  any  integer  n,  i.e.,  the  input  plane  field  is  replicated  at  each 
of  the  Fourier  image  planes  z  =  nz^.  This  phenomenon  was  first  investi¬ 
gated  theoretically  for  unbounded  periodic  gratings,  and  the  imaging  dis¬ 
tance  was  first  given  by  Rayleigh  [14],  Later  treatments  have  been  given 
in  [13],  [15]-[17],  Because  of  the  mathematical  equivalence  of  a  wave¬ 
guide  with  a  periodic  system,  this  imaging  also  occurs  in  waveguides--a 
fact  apparently  first  noticed  by  Rivlin  and  Shul'dyaev  [18]  and  discussed 
at  length  by  a  number  of  authors  [10],[19]-[23], 

An  even  more  interesting  occurrence  shows  up  at  z  =  Zp^,  where 


z 


pq 


(21) 


and  p  and  q  are  some  positive  integers.  Let  us  consider  the  sum 

2 

°°  -iimrx/a  +  2-rcim  z  /z, , 

V*>  -J'  pq 

m—  cct 


(22) 
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Let  m  =  p£  +  r,  where  r  runs  from  0  to  p-1,  and  express  (22)  as  a 
double  sum: 

Q  (x)  =  l  e-i™(p£+r)/a  +  2Triq[p£2+2£r  +  r2/p] 

pq  £=-«  r=o 

OQ.  V 

P£  g-t^rx/a  +  2irtqrZ/p  (  l  e-VFTp£x/aJ 

r=Q  £  — -  aa 


(23) 


However,  the  summation  in  parenthesis  is  nothing  more  than  Q-j  1  ( px ) »  although 
now  the  argument  can  range  not  just  from  -a  to  +a  (or  0  to  2a),  but  from 
-pa  to  +pa  (or  0  to  2pa).  Making  use  of  the  formula  [24]: 

l  e2ir1mx/d  =  d  l  fi(x-nd)  (24) 

m*  ~  »  n=  -  oo 


for  any  positive  d, 

P-1 

V*1  ■  T  Jo  *' 


iTrrx/a  +  2iriqr  /p 


l  Six 


2na, 


n=  -00 


(25) 


From  (18),  then,  we  have 
-ikz_ 


G0(x’X';Zpq)=e 


pq 


l  c  (p.qKSU-x'-^-fiU  +  x'-^-)] 

n=  v  H 


2na 


(26) 


where  the  coefficients  cp  are  given  by 


(p,q)  -  1  V  e2’Ii'r(r<»+n)/P 
1  p  r=0 


(27) 


For  p  =  1,  only  one  of  the  delta-functions,  6(x-x‘),  is  nonzero 
in  the  range  0_<(x,x’)  _<  a,  and  we  recover  the  single  Fourier  image 
described  before.  If  p>l  on  the  other  hand,  more  of  the  delta-functions 
in  (26)  may  appear  in  this  range.  Each  one  contributes  to  (13)  a  replica 
of  the  input  field  which  is  shifted  by  some  amount  in  the  x-direction. 
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and  whose  amplitude  is  |cn(p,q)|  times  that  of  the  original  image.  Any 
terms  arising  from  the  terms  5,(x  +  x'  -2na/p)  are  inverted  as  well.  Images 
of  this  type  have  been  called  Fresnel  images. 

For  an  input  function  not  symmetric  with  respect  to  the  center  of 
the  guide  x  =  a/2,  we  have  depicted  the  various  images  along  with  their 
(complex)  amplitudes  in  Fig.  2  for  z^ ,  z31  and  z^,  The  phase  factor 
exp(-ikz)  is  omitted  in  this  figure.  Pigure  2(a)  illustrates  the  input 
function  Ey(x,0).  At  z = z^  ,  there  is  only  one  image,  which  is  inverted 
with  respect  to  the  original.  But  is  also  a  further  180°  out  of  phase  with 
respect  to  the  input  after  the  factor  exp(-ikz)  is  accounted  for  (Fig.  2(b)). 
At  z  =  z.^  (Fig.  2(c)),  the  situation  is  more  involved.  One  of  the  images 
is  an  unchanged  replica,  reduced  in  amplitude  by  1//3"  and  phase  shifted 
by  90°.  The  other  two  image  terms  have  "broken  up"  and  rearranged  the 
original  pattern.  Between  them,  they  both  contain  two  complete  replicas 
of  the  original,  and  again  the  amplitude  of  each  is  reduced  by  1//3? 

Note  that  this  amplitude  reduction  is  consistent  with  the  fact  that  all 
of  the  power  of  the  original  pattern  must  be  divided  between  the  three  split 
images.  At  z  =  z^,  there  is  both  an  erect  and  an  inverted  image  (Fig.  2(d)). 

For  the  special  case  of  a  symmetric  excitation,  Ey(x,0)  =  Ey(a-x,0), 
z^,  z^  and  indeed  Zg^  all  reproduce  the  original  input  function  (Fig. 3(a)). 
If  the  input  function  is  a  beam  of  sufficiently  narrow  width,  we  can  recog¬ 
nize  three  essentially  distinct  images  at  z^  (Fig.  3(b)).  At  other  image 
planes,  similar  conclusions  hold. 

The  Fourier  and  Fresnel  images  allow  us,  in  principle,  to  compute 
the  field  in  the  waveguide  at  any  point  z  (and  any  arbitrary  value  of 

z  can  be  approached  as  nearly  as  desired  by  such  a  point).  This 
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procedure  could,  however,  require  a  large  value  of  p,  and  hence  an  inordi¬ 
nately  large  number  of  image  terms,  resulting  in  a  method  which  is  no 
more  efficient  than  the  modal  approach.  In  the  next  section,  we  will 
show  how,  for  a  certain  specific  type  of  excitation,  efficient  field  compu¬ 
tation  can  be  carried  out  for  any  value  of  z,  using  only  a  relatively 
small  number  of  images. 

Let  us  emphasize  in  closing  this  section  that  the  imaging  phenomenon 
results  from  the  collective  interference  of  the  mode  sum  (14),  or  alter¬ 
natively  from  the  interference  of  the  series  of  MraysM  (plane  waves) 
represented  by  (18),  depending  on  one's  preferred  physical  picture.  The 
only  assumption  involved  is  the  paraxial  approximation,  and  in  the  case 
of  the  parallel-plate  waveguide,  no  approximation  of  the  mode  functions 
themselves  is  needed.  Tn  the  case  of  a  dielectric  waveguide,  some  small 
higher-order  corrections  to  the  mode  functions  will  be  needed  (see 
Appendix  A), 
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IV.  PROPAGATION  OF  A  GAUSSIAN  BEAM 


Consider  the  initial  field  distribution 

-(x-x  )2/2w2 

Ey(x,0)  =  e  0  0  (2 

i.e,,  a  Gaussian  beam  centered  at  x.  with  waist  parameter  w  . 

0  0 

Let  0<xQ<a,  and  assume  that  the  "tails"  of  the  beam  are  negligible  at 
the  walls  of  the  guide: 


wo  ^  xo  ; 


w„  <<  a  - 
o  o 


In  addition,  we  also  suppose  that  the  beam  is  well  collimated,  kwQ  » 1 
Under  these  conditions  and  the  paraxial  approximation,  (13)  can  be 


replaced  by 


r  -(x'-xjW 

r(x,z)  =  e  0  G0(x,x';z)dx 


We  wish  to  evaluate  (29)  for  arbitrary  (not  necessarily  rational)  va’ues 

Of  Z/Z-,  y 

To  do  this  it  is  convenient  to  express  GQ  in  terms  of  the  Jacobian 
theta-function  defined  by  Whittaker  and  Watson  [25]  as: 

<?3(z|T)  =  l  em27riT  +  2miz  (30) 

m=  -oo 

(the  argument  z  is  conventionally  used  in  this  connection  and  should  not 
be  confused  with  the  cartesian  coordinate  z  used  above).  This  function 
has  a  wealth  of  useful  properties  which  we  summarize  in  Appendix  B.  The 
theta-function  appears  in  solutions  of  other  parabolic  equations,  such  as 
the  heat  equation  [28].  In  the  paraxial  parabolic  equation,  however,  the 
arguments  z  and  x  are  both  real  (see  below)  and  0^  must  be  treated  as 
a  generalized  function. 
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Using  (30),  GQ  becomes 


r  («  v'.,\  -  1  -tkz  I 
c>  "  2d  ^  ^ 


2z 

Z11 


(31) 


From  (29),  then 

,  >  e'1kz  r  -Cx'-xo>2/Zw" 

Ey(x>zl  =  ST-  '  e 


With  the  help  of  (B.ll)  we  obtain 


°  |  it,)  - 


it 

(32) 


Ey(x,2)=fyfe-" 


L(’,(x'x°> 

2z 

,  2 

+  l7TWo  \  J 

2z  + 

'""I  (l 

pi  2a 

Z11 

1 Ti 

3l  2a 

zn 

2a2  J 

(33) 


Let  us  first  consider  the  focussing  relative  to  z  Let 
z  =  qz^  +  Az,  where  -zp2  <  AZ  <z^/Z  and  q  is  an  integer.  From 
(B.3)  and  (B.8),  we  have 


Ux,z)  =7 


w 


fftiy 


o  -ikz 


<0)2/2f2( 


('  a(x-xQ) 

2ia2  ) 

'  f2(Az) 

irf2(Az)  ' 

(34) 


1 2  # o^2 / 


p(x+xQ) 

2ia2  )j 

'  f2(Az) 

irf2( Az)  '] 

where  we  have  defined  a  "complex  waist  parameter"  f(Az)  as 

.2 


’(Az)  =  w2[] 


1  zt;]  =  wo  ■  1'Az/k 


4a‘ 


ttw„  IT 
o 


(35) 


Using  (B.12),  we  can  further  reduce  (34)  to 


w  _ikz  »  f-(x-xQ+ 2ma)2/2f2(Az)  -(x+xQ+2ma)2/2f2(Az)l 

’(x’z)’fl4lTe  „.Me  'e  J 

k  (36) 


Now 
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_ 1_ 

2f2(Az) 


2w2(Az) 


iAz/kw2 

2w2(az) 


where  the  waist  size  w(a z)  is  given  by 


(37) 


w2(az)  *  w2  +  (Az/kwQ)2  (38) 

so  we  can  see  that  (36)  represents  an  infinite  series  of  Gaussian  beams, 
each  broadened  from  its  focal  plane  z  = qz^  as  if  it  were  propagating  in 
free  space.  We  illustrate  this  in  Fig.  4.  As  az  increases,  more  and 
more  of  the  "image"  beams  contribute  significantly  to  the  field  in 
0<x<a.  At  az  =  zn/2  (say),  the  waist  size  has  become 

w2(Zll/2)  =  w2  +  (2a2/m*0)2  *  4a4/n2w2  »  a2 


because  our  assumptions  about  the  beam  imply  a  »  wQ.  Thus  we  may 
require  quite  a  few  terms  of  the  image  series  (36)  in  order  to  compute 
the  fields  at  certain  values  of  z. 

For  such  values  of  Az,  we  seek  to  improve  the  efficiency  of  our 
scheme  by  transforming  (33)  in  a  somewhat  different  manner.  Applying 
relation  (B.9),  we  find  that 


Ey(*.z)  -  ^/fVikz  £  e’1r2ZZ/Z"  -,2r2"2/2a2 

^  *  ir  l4]|  p2  £  ♦  ip2 

3  za  zn  2gz  i  zn 

.  1r  lM]|p2|z  i 

U  4a  zn  2a4  I  zn 


(39) 
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for  any  specified  positive  tnteger  p.  Now,  let  z  =  z  +  Az  ,  where 

rM  r 

now  -z^/2p  £  AZp  <  z^/2p.  The  periodicity  properties  (B.3)  and  (B.4) 
allow  us  to  replace  z  by  AZp  in  the  arguments  of  the  theta-functions 
in  (39).  Subsequently  applying  (B.8)  results  in 


Ey<‘-2)-OT%Te'ik2  t* 


2uir  q/p 


r=0 


-(x-xQ)2/2f2(Azp)  ,  f  ia(x-xft) 


_ o  _  mr 

pf2(Azp)  P 


2ia2 


irp2f2(Azp) 


)- 


-(x+xQ)2/2f2(Azn)  _  / ia(x+xj 


'M 


_ ITT 

3V2(azp)  p 


irp 


■t2  4J./J 


P 

(40) 


where  f  is  defined  in  (35)  as  before.  For  p  =  l,  this  reduces  to  (34) 
as  expected.  Finally,  through  the  use  of  (B.12)  we  obtain 


c  lv  _ _ o  -ikz  1 

V  ’  "f(AZp)  e 


L^1  e2^r2q/pJ  J  e-2,in,r/Pre-(x-V2ma/P>2/2f2^p> 

0  m=  -oo  i 

- ( x+xQ+2ma/p ) 2/2f 2 ( Azp  )J 


(41) 


Clearly,  equation  (41)  represents  a  string  of  Fresnel  images  at  zp^, 
broadened  by  their  additional  propagation  distance  Azp  as  evidenced  by 
the  factor  f(Azp) (compare  (26)-(27)).  A  slight  rearrangement  of  (41) 


yields  a  single  summation  (letting  m-*--n): 


Ey(x>2) ‘fT^T  e"k2  jLcn(p-,,)^ 


-(x-xQ-2na/p)2/2f2(Azp) 


-e 


-(x+x0-2na/p)2/2f2(Azp)| 


(42) 


where  cn(p,q)  is  given  by  (27). 

How,  then,  is  p  to  be  chosen?  Our  goal  is  to  minimize  the 
number  of  terms  of  the  series  (42)  required  to  give  a  specified  accuracy. 


If  we  elect  to  truncate  the  series  at  n =  ±N  such  that  the  terms  at  this 

-ft 

point  have  magnitude  less  than,  say,  e  ,  N  will  be  the  next  larger 
integer  than 

2pw(Azp)/a 

where  w  is  given  by  (38)  and  we  have  examined  the  "worst  case"  when 
x  =  xQ  in  the  first  term  of  the  summand  in  (42).  As  Azp  varies  from  0 
to  z^/2p,  this  quantity  varies  from  2pwn/a  to 


o 

2  2  4 


-  —  (l  +  — ~T—  ) 
7T  wQ  V  4a4  } 


If  p  <  a/2wQ,  N  will  vary  between  a  minimum  of  1  and  a  maximum  of 

2  2 

about  4a/Trwo.  If,  on  the  other  hand,  p  >  2a  /ttwq,  then  N  is  at  least 
4a/TTwQ  as  a  minimum,  and  at  least  4/*^ a/TTWQ  as  a  maximum.  Clearly 
then,  such  large  values  of  p  are  undesirable  for  efficient  field 
computations  (because  more  images  of  the  Gaussian  beam  are  fitted  inside 
the  waveguide  than  can  be  accommodated  without  severe  overlap),  and  there 
is  no  apparent  reason  to  let  p  exceed  a  value  on  the  order  of  a/wQ. 

On  the  other  hand,  for  a  specific  value  of  z,  Azp  may  be  much  closer 
to  zero  for  a  smaller  value  of  p,  and  in  such  cases  the  smaller  choice 
of  p  might  be  desirable.  However,  a  single  choice  of  p  for  all 
values  of  z  results  in  a  simpler  computer  program,  and  it  is  seen  from 
the  foregoing  discussion  that  a  value  of  about  a/wQ  is  nearly  optimal. 
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V.  NUMERICAL  RESULTS 


Numerical  results  for  Q^z £Zg1 


were  computed  for  a  symmetrical 


Gausstan  beam 


Ey(x,0)  =  e 


-(x-a/2)2/2w2 


for  a  waveguide  with  ka  =  973.4,  kwQ  -  11 Q. 3.  Computations  were  made 
using  both  an  exact  mode  series  (cf.  (14))  and  the  hybrid-image  represen¬ 
tation  (42),  The  two  methods  gave  graphically  indistinguishable  results 
over  this  range,  which  are  displayed  in  Fig.  5.  Note  that  5  or  6  is 
the  maximum  distinguishable  number  of  beams  in  this  case,  and  is  a 
suitable  choice  for  p. 

Figure  6  shows  that  for  very  large  values  of  z,  the  accuracy  of 
the  paraxial  approximation  has  begun  to  deteriorate  slightly.  In  a 
future  paper,  we  will  examine  how  to  obtain  closed-form  expressions  for 
this  correction,  but  meanwhile  we  note  that  even  for  z  =  250z^  ,  the 
accuracy  of  the  paraxial  expression  is  quite  good.  If  we  take  this 
example  to  model  an  optical  waveguide  with  a  =  100y,  wo  =  7.07p,  the 
distance  250z^  represents  about  a  50m length  of  waveguide  which  of 
course  is  a  huge  number  of  wavelengths.  We  thus  see  that  the  paraxial 
approximation  is  capable  of  excellent  accuracy  over  modest  lengths  of 
waveguide,  even  at  optical  frequencies. 
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VI.  CONCLUSION 

Ke  have  described  a  hybrid  technique  for  computing  the  fields  of 
a  paraxial  beam  propagating  in  a  multimode  waveguide  for  very  long 
distances.  The  method  relies  on  the  periodic  Fourier  and  Fresnel 
imaging  properties  of  the  guide,  and  is  highly  efficient  for  beams  of 
moderate  width  compared  to  either  full  modal  or  ray  approaches.  Numer¬ 
ical  comparisons  have  confirmed  the  accuracy  of  this  method. 

This  approach  should  be  susceptible  to  generalizations  in  several 
directions.  For  a  beam  with  oblique  incidence  (as  in  [6]),  the  parabolic 
approximation  of  section  II  can  be  modified  when  the  dominant  propagation 
factor  in  the  z-dfrection  is  other  than  exp(-ikz)  as  indicated  in  [12]. 
Higher-order  corrections  to  the  paraxial  solution  for  a  dielectric  slab 
similar  to  those  for  the  metallic  guide  given  in  Section  II  can  also  be 
obtained.  Slowly-varying  guide  widths  and  inhomogeneous  refractive 
index  profiles  should  also  be  tractable  by  similar  methods.  Investigations 
into  these  areas  are  currently  being  made,  as  is  the  generalization  to 
waveguides  of  circular  symmetry. 
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APPENDIX  A 

THE  PARAXIAL  APPROXIMATION  FOR  A  DIELECTRIC  SLAB 


Consider  the  step-index  dielectric  slab  waveguide  shown  in  Fig.  A.l. 


The  slab  has  thickness  b  and  refractive  index  nQ.  The  cladding  index  is 
n^,  and  both  media  are  assumed  to  be  nonmagnetic.  The  TE  modes  for  this 
waveguide  have  the  field  distribution  [1] 


e-1kaz  ft  e 


Iz  2 

•/a  -n. 


x  <0 


E  = 

y 


sin  (j)  (a) 

1'*<aZ  A  sin  jkv^  - x  + 

•ikaz  f  /?  2  1  -k^-n^x-b) 

lkaZ  A  sinjk/r  -a  b  +  <j>(a]|  e  1  x  >1 


(A.l) 


0  <  x  <  b 


where  A  is  an  arbitrary  constant  amplitude.  Here  k  is  the  wavenumber 
of  free  space,  ka  is  the  propagation  constant  of  the  mode,  and  (p(a ) 
is  the  phase  shift  associated  with  the  Goos-Heinchen  effect: 


where  V 


4>(a)  =  sin’^kb  /n^ -ot2/vJ 

is  the  so-called  normalized  frequency  : 


(A. 2) 


(A. 3) 


Note  that  V  >>  1  for  a  highly  multimode  guide. 

The  characteristic  equation  which  determines  the  eigenvalues  a 
obtained  by  requiring  Hz  to  be  continuous  at  x  =  b: 

sin  jk/i^  -  a2  b  +  2<j>(a)l  -  0 


is 


(A. 4) 


The  paraxial  approximation  to  these  modes  (o  =  nQ)  is  found  as  in  [10] 


A-2 


by  reckoning  4>(a)  to  be  small.  From  (A. 2)  and  (A. 4),  we  then  obtain 
approximately 


sin 


=  0 


or 


a 


m2  2 
m  it 


k2b2(l  +  2/ V ) 2 


(A. 5) 


i.e.,  the  propagation  constants  for  a  parallel-plate  waveguide  of 
slightly  larger  width  a=b(l  +2/V).  The  corresponding  field  within  the 
slab  is,  from  (A.l ), 


-ika  z 
m 


A  sin 


Lk^  -  a2  (x  +  b/V) 


0 


0  <  x  <  b 


Here  we  have 


oc  =  n 
m  o 


2  2 
m  ir 


2nQk2b2(l  +  2/V ) 2 


+  0 


44 

m  it 


k2b2V3(l  +  2/V)5 


4  4 
m  it 


L  k4b4(l  +2/V)4" 


(A. 6) 


(A. 7) 


For  most  optical  waveguides  the  first  error  term  is  likely  to  be  the 
larger,  but  in  any  case  we  require 


k2b2V3  »  1 


4  4 

and  k  dh  »  1 
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APPENDIX  B 

PROPERTIES  OF  x|(z|t) 

In  this  Appendix,  we  derive  a  number  of  useful  properties  of  the 
theta-function  i?3(z[x),  defined  by 

>1(z|t)  =  l  em2lliT  +  2miz  (B.  1 ) 

J  m  =  -  00 


where,  for  our  purposes,  z  is  an  arbitrary  complex  number,  whi^e  x 
lies  in  the  upper  half-plane  Im(x)  >0.  If  x  lies  on  the  real  axis, 
we  likewise  require  that  z  be  real,  and  in  this  case  i^(z|x)  must  be 
treated  as  a  generalized  function. 

From  Whittaker  and  Watson  [25],  we  can  obtain  a  number  of  periodicity 
and  parity  relations: 

^3(z|t)  =  ^3(-z|x)  (B.2) 


■^3(z|x  +2n)  =  ^3(z|t) 

xSj(z  +  nir[x)  =  z  j  x ) 

$j(z  +mrx|x)  =  e-7r1n^T  "  2inZi^(z|x)^ 

All  four  relations  are  easy  consequences  of  the  definition  (B.l). 
also  interesting  to  note  that  satisfies  the  parabolic  equation 


l  n  =  0 ,  ±1 ,  ±2 , . . 


(B.3) 

(B.4) 

(B.5) 


It  is 


XT  i 

4 


a^UI-r) 


(B.6) 
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which  is  obtained  from  the  Helmholtz  equation  in  the  paraxial  approximation, 
and  is  similarly  easily  verified  (cf,  eqn,(8)). 

Per  real  z,  at  t»0  (as  a  generalized  function)  can  be  evalu¬ 
ated  as  [24], 126]: 

^3(z|0)  =  l  e2miz  =  tt  l  6(z-mr)  (B.7) 

m=  n=  -°° 


Another  useful  relation,  which  holds  for  general 
obtained  from  Jacobi's  imaginary  transformation  [21]: 


z  and  t,  is 


(B.8) 


This  relation  can  be  verified  using  the  Poisson  summation  formula  [24]. 

The  identities  which  form  the  basis  for  the  "image-splitting" 
properties  of  the  theta  -function  can  be  deduced  from  a  more  general 
expression  given  by  Krazer  [27].  These  relations,  which  might  be 
referred  to  as  modular  relations,  are 


lA(z|x)  -  e7rir2'r+2irZ  1$,(p[z+irn:] |p2i)  (B.9) 

J  r=0  J 


Equation  (B.9)  is  verified  by  writing  m  =  pJt+r,  r  =  0,1 , . . . ,  (p-1 )  in  (B.l); 
(B.10)  follows  by  substituting  (B.l)  into  the  right  hand  side.  Actually 
(B.9)  and  (B.10)  can  also  be  derifed  from  each  other  using  (B.8)  as  well. 
An  integral  which  the  authors  could  not  find  in  the  literature 


is  given  below: 


B-3 


“  e-ax2tbx  v«(x|T)dx  ,  i$(&|t+£  >  (B.H) 

-  oa 

The  derivation  is  straightforward,  proceeding  by  integrating  (B.l)  term-by- 
term. 

Finally,  we  note  that 

li(1z|lT)  -  e2^/lTT  J  e-(z+mxT)2/xt  (8.12) 

m=  -  00 

i.e.,  can  be  related  to  a  string  of  displaced  Gaussian  functions. 


Fig.  1:  Parallel-plate  waveguide. 

Pig.  2:  Imaging  of  a  non symmetrical  field  distribution. 

(a)  Input  function  Ey(x,0). 

(b)  Ey(x,z21)  =  Ey(a  -  x,0)  (Inverted  image). 


(c)  The  three  components  of  Ey(x,z^); 

i  e-W6 

Ey(x’z3i}  =  ^Ey(x,0)  +  ~7T 


-Ey(^-  -  x,0)  ; 
Ey^x  "  T  ,0^ 


■  i  7r/6 


+  fi- 


Ey(x  +  °)i  CUx^ 


j-E  (y  -x,0)  ;  x^a 

(d)  The  two  components  of  E  (x,z^); 

Ey(x,z41)  -  Ce1  71/4  Ey ( x , 0 )  +  e3rri/4Ey(a  -  x,0)]/  . 


0<  X-S  ~ 

2a 

■y«  x<  av 


Fig.  3:  Imaging  of  a  symmetric  field  distribution. 

(3)  |Ey(x,0)|  =  |Ey(x,z2i)|  =  |Ey(x,z^.|)|  =  |Ey(x,Zg^)|  . 

(b)  The  three  components  of  Ey(x,z^). 


Fig.  4:  Overlapping  of  original  and  image  beams  in  overmoded  waveguide. 

Fig.  5:  Evaluation  of  power  distribution  of  Gaussian  beam  over  1/8  -cycle: 

ka  =  973.39;  kwQ  =  110.3.  (a)  z/z^-  0;  (b)  3/512;  (c)  8/512; 

(d)  13/512;  (e)  16/512;  (f)  22/512;  (g)  28/512;  (h)  32/512; 

(i)  55/512;  (j)  64/512. 

Fig.  6:  Exact  and  approximate  power  patterns  at  large  axial  distance; 
ka  =  1538.23;  kwQ  =  174.3;  z/z^  =  250.119. 


Fig.  A.l:  Dielectric  slab  waveguide 


